Direct reconstruction of dark energy 
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An important issue in cosmology is reconstructing the effective dark energy equation of state 
directly from observations. With so few physically motivated models, future dark energy studies 
cannot only be based on constraining a dark energy parameter space. We present a new non- 
parametric method which can accurately reconstruct a wide variety of dark energy behaviour with 
no prior assumptions about it. It is simple, quick and relatively accurate, and involves no expensive 
explorations of parameter space. The technique uses principal component analysis and a combination 
of information criteria to identify real features in the data, and tailors the fitting functions to pick 
up trends and smooth over noise. We find that we can constrain a large variety of w(z) models to 
within 10-20% at redshifts z < 1 using just SNAP-quality data. 



Introduction The dark energy crisis in cosmology 
highlights our incomprehension of what the universe ac- 
tually consists of. Usually characterised by an effective 
equation of state function w in which we hide our lack of 
understanding, an important goal over the coming years 
will be to try to understand this variable as a function 
of redshift, giving w(z). Theoretically we have little to 
go on, other than the cosmological constant which has 
w = — 1 for all time. Alternatives such as quintessence 
or modified gravity theories make predictions about how 
w diverges from —1 yet are often forced to parameterise 
free functions [1 j. More radically, proposals which modify 
the radial distribution of matter on Hubble scales provide 
no a priori constrains on what the effective equation of 
state (defined by matching up the distance indicator with 
an FLRW model) could be (2J. So, the forward problem 
of parameterising models, matching to data and discard- 
ing the fits which are poor or over-parameterised, suffers 
from a profound arbitrariness when we try to interpret 
the errors: using simple smoothly varying functions of z 
severely limit the departures from cosmological constant 
to a small range of models, but if we add more freedom to 
w the errors grow uncontrollably. Without encapsulating 
the behaviour of w(z) which may actually exist, what are 
we really constraining? Would a 'backwards' method be 
better? Can we instead reconstruct w from observations 
directly? 

The dark energy equation of state is typically 
(re) constructed using distance measurements as a func- 
tion of redshift. The luminosity distance may be written 

as d L (z) = Hol/~-n k sin (y-^fc Jo w[p)) > where H(z) 
is given by the Friedmann equation, H(z) 2 = Hq {Q m (l-\- 

z f + n k (i + zf + (i - n m - n k ) ex P [3 ; 2 

where Hq = H(0) and fL m ,k are the normalised den- 
sity parameters. A common procedure is to postulate 
a several parameter form for w(z) and calculate (Il(z). 
The most promising of these approaches uses a princi- 



pal component analysis to construct the 'optimal' ba- 
sis functions for w(z) based on the data [3]. An alter- 
native method is to reconstruct w(z) by directly recon- 
structing the luminosity-distance curve. Writing D(z) = 
(Hq/c)(1 + z)~ 1 di J (z) as the normalised comoving dis- 
tance, we have [5HZ]- 

w(z) = {2(1 + z)(l + n k D 2 )D" - [(1 + z) 2 n k D' 2 
+2(1 + z)Q k DD' - 3(1 + Sl k D 2 )]D'}/ (1) 
{3{(1 + z) 2 [Q k + (1 + z)n m ]D' 2 - (1 + n k D 2 )}D f }. 

Thus, given a distance-redshift curve D{z), we can re- 
construct the dark energy equation of state, assuming 
we know the density parameters Q m and Q k . Different 
methods for doing this involve smoothing the data to give 
D(z), or parameterising D(z) by a function; see [4 for a 
comprehensive review, and [9HT5] for alternative model 
independent approaches. 

The direct reconstruction method is unstable because 
of the two derivatives of the observed function in eq. 
|2j requiring the fitting function to accurately capture 
the slope and concavity of luminosity distance curve. 
This means that differences between the true underlying 
model and the fitted function due to the choice of param- 
eterisation, are amplified drastically when reconstructing 
w. Furthermore, w is constructed from a quotient of 
functions which need to balance to obtain the correct w. 
However, there must exist a set of 'optimal' basis func- 
tions with which achieving this balance becomes possible 
and direct reconstruction feasible; in essence this is the 
same as specifying 'the correct' parameterisation for dark 
energy. Can we find such a basis? 

At first sight it seems not: If we have no inherent intu- 
ition of the true w, how can we possibly guess the right 
form for D(z)7 Of course, a large polynomial expansion 
would work but at the expense of ludicrous errors. Fur- 
thermore, one can easily achieve a fit that is too good: 
a x 2 i s l ess than the value the actual underlying model 
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FIG. 1: Extracting w = — 1. The first 10 eigenf unctions for A are shown left. We take two separate reconstructions, using 
M — 4 and M — 5 (N = 10), and generate the errors on /i res id via a monte carlo (with 100 samples each for clarity). The 
two reconstructions are combined with equal weight to give the 1-a errors. The choice of number of eigenfunctions is not yet 
favourable; beyond z ~ 1.2 (not shown) the reconstruction is hopeless for these eigenfunctions. 



would produce. Over-fitting to noisy data translates into 
wild behaviour in w. 

Here we present a method to find well adapted ba- 
sis functions for fitting D(z) and a simple way to con- 
struct the errors on w. We first calculate the residuals of 
the measured apparent magnitude around some fiducial 
model. Assuming some basis functions, we then calculate 
the principal components of the residuals, and use those 
fixed principal components to provide a measure of the 
true D(z). Given the many possibilities, a combination 
of information criteria are used to select which encom- 
pass the information present in the data. Folding the 
errors together appropriately produces a non-parametric 
method which reproduces w(z) together with an effec- 
tive '1-cr' confidence measure. A non-parametric method 
doesn't produce confidence limits; rather a confidence in- 
terval of x% must be expected to trap the correct value 
x% of the time [16]. 

Basis functions for distances Let's assume we have 
Nd data points distributed at Z{ for the distance modu- 
lus fx = 5 log 10 oIl(z) + 25, with Gaussian errors which 
are independent. We create p resid = /i data - /ifiduciai 
where the fiducial model is some predefined model, such 
as flat LCDM, an empty model, or EdS. A good choice is 
the best-fit LCDM model, being consistent with current 
data. Our goal is to construct /i res id and two deriva- 
tives as accurately as possible. We choose a set of pri- 
mary basis functions p n (z) such as p n (z) = z n ~ x and use 
the function Y^i=\ a nPn(z) as the basis to fit /i re sid(^) 
to data, for a fixed N . This is a linear fit, so can be 
done easily, and the covariance matrix C calculated al- 
gorithmically without having to explore a complicated 
parameter space, or make assumptions about Gaussian 



errors on the parameters. Now we perform a princi- 
pal component analysis on the fit to find the best ba- 
sis of functions as follows. We diagonalise the inverse 
covariance matrix, and create a matrix of eigenvectors 
E. We then order them according to decreasing eigen- 
value; i.e., so that the n'th column of E, e n , is the eigen- 
vector with the n'th largest eigenvalue, etc. Write e mn 
as the m'th component of the n'th eigenvector. Define 
the eigenfunctions P n (z) = Y,m=i e mnPm(z), for n E 
[1,7V] which now form a new family of basis functions 
which are suitably adapted to the data. In matrix 
form P = E T p <^> p = EP. These functions are 
now orthogonal with respect to the errors on the data: 
^2f d a~ 2 P m (zi)P n (zi) = if m ^ n. In other words, the 
Fisher matrix in this basis is diagonal. We can normalise 

the eigenfunctions P n (z) = P n {?) I '\jYn d a i~ 2 p n(zi) 2 , m 
which case the Fisher matrix will be the identity ma- 
trix. If we now refit to the data with new parameters for 
these basis functions, the covariance matrix is also (very 
nearly) the identity matrix. For a fixed iV, we may be 
interested in the first M < N eigenfunctions which en- 
compass the dominant features in the data, and throw 
away the higher ones which contain noise-induced oscil- 
lations. 

Practicalities and a test case - can we recover A? To 
test the method, we construct a hypothetical data set in 
line with expectations from the SNAP experiment, con- 
sisting of 2000 type la supernovae (SNIa) measurements 
evenly distributed in the redshift range z = 0.08 — 1.7, 
and 300 local SNIa in 0.03-0.08 [T7]. The statistical un- 
certainty is conservatively estimated as o- mag = 0.15mag 
and we include systematic errors, as a linear drift from 
o~ S ys = — 0.02 [17]. We assume a true underly- 
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FIG. 2: Reconstruction of w — — 1 using the CIC criteria, Eq. Q, with s — 0.2. All [TV, M] values satisfying Eq. Q are monte 
carlo-ed using 200 runs each (left); the results are bundled into a single probability distribution (middle) from which we infer 
1-cr errors at each redshift, shown left. For this example, the CIC selects for M — 2,3 for each TV. We show two probability 
distributions for TV = 10 (middle), as well as the combination for all TV, M used. It is at the stage of bundling up the different 
TV, M eigenmodes into one probability distribution where the method becomes non-parametric. The combination of different 
eigenmodes leads to slight non-Gaussianity of the distribution. In this example, s can be increased all the way up to 1 which 
shrinks the errors even further, and allows tight constraints on w over the full range of z (right). 



ing w and some parameters, £ tme = {^ m , fife, h} true = 
{0.3,0.0,0.65}. Our goal is to recover w true (z) in a ro- 
bust way. We assume £ fiducial = £true f or now ^ anc [ ^ a ]^ e 

^fiducial = — 1- The actual model used for //fiducial is 
re-incorporated when evaluating w(z), but an incorrect 
choice of the parameters for £ fiducial are subsumed by the 
usual uncertainties on those parameters. If we pick a 
fiducial model which is reasonably close to the true un- 
derlying cosmology, the residuals are predominantly noise 
and variations in w, requiring fewer fitting parameters 
giving typically smaller errors. Errors from the incorrect 
choice of £ will be considered elsewhere but are standard. 
Here, we consider only the reconstruction errors for clar- 
ity. 

We used a variety of primary basis functions such 
asp n (z) e {z n -\[z/(l + z)] n - 1 ,[l/(l + z)] n - 1 }, with 
similar results, though we achieve smoother reconstruc- 
tions as we move left to right in this list; we present 
our results using the middle one. For a fixed TV the 
eigenfunctions range from smooth with no turning points, 
to very oscillatory; typically, the n'th eigenmode crosses 
zero n — 1 times. Roundoff error can cause problems for 
TV > 10, which is signalled by C having non-zero diago- 
nal elements for the normalised basis. (Strictly speaking, 
C differs from the identity matrix slightly, since we use 
one realisation of the data: we would expect a change of 
roughly l/y/N^ for the diagonal elements.) For TV > 10 a 
singular valued decomposition could be used for the fits. 
In Fig. [T] we show the first 10 normalised eigenfunctions 
when the underlying model is w — —1, with TV = 10. 



The first M of these eigenmodes are used as the new ba- 
sis functions for /i re sid, and these are fitted to the data. 
In the normalised basis the errors on the parameters are 
all Gaussian and unity (up to ~ l/y/Nd). Since the pa- 
rameters are uncorrelated, the errors may be propagated 
into /i r esid by a simple monte carlo for each parameter. 
These curves then give a family of w (z) curves from which 
the 1-cr error may be given. Errors on £ may be folded 
in at this stage and lead to larger error bars, though we 
don't investigate this here. We show an example of this 
reconstruction procedure in Fig. [I] using M = 4 and 5, 
mixed together to form one set of error bars. There is 
nothing to say that these are good choices of M, an issue 
we explore now. 

Selection criteria The number of eigenfunctions to 
use in the final reconstruction is critical as it deter- 
mines the accuracy and size of the errors bars. Consider 
a subset of [TV, M] with TV £ [2,10] and M e [2,7V]. 
Each choice [TV, M] will give a particular reconstruction 
of w(z), together with some errors. In the case where 
M — TV, the original error bars are recovered and no in- 
formation is thrown away. Reducing M is accompanied 
by a reduction in the error, but an increased chance of 
getting w(z) wrong. How do we select the 'correct' set of 
eigenfunctions? Choosing the combination [TV, M\ lead- 
ing to the smallest x 2 runs the risk of overfitting to noise. 
The Risk may be used [3j [16] , but requires the knowledge 
of the underlying value of w(z) a priori. 

Instead, we use a mixture of Akaike and Bayesian in- 
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FIG. 3: Reconstruction of evolving dark energy. With a choice of s = 0.2 and n = 5 we can reconstruct w to z ~ 1 with 
~ 10-20% accuracy. 



formation criteria [18]: 

AIC = Xmin + 2M, BIC = x2u„ + MlnJV d , (2) 

where smaller values are assumed to imply a more 
favoured model. The utility of these criteria over the 
Risk is that they are computed without knowing the un- 
derlying solution. We evaluated AIC and BIC values 
corresponding to each [TV, M\ combination for a number 
of test cases and found that minimizing these two criteria 
lead to dramatically different reconstructions that were 
usually not optimal. Typically, models with a very low 
BIC are too smooth with tight error bars, while those 
with low AIC values are often too oscillatory and have 
large errors. A more adaptable requirement uses a com- 
bined information criteria which we define as: 

CIC = (l-s)AIC + sBIC, (3) 

where s takes us from conservative models when s = 1 to 
more wild models when 5 = 0. The parameter s thus me- 
diates the competition between a smoother reconstruc- 
tion (in which BIC is minimized) and one that is more 
featured (in which the AIC criterion is smallest). 

But there is no reason to select one particular recon- 
struction; the minimum CIC is still no silver bullet. We 
find a successful strategy is to select different [TV, M] 
choices which are near the best values of CIC, for a given 
5, and amalgamate them at the monte carlo stage when 
we compute the errors. We weight each [TV, M] choice 



equally. In this way, we reduce inherent bias which ex- 
ists in any particular choice of [TV, M], even after the 
principle component analysis. 

After experimenting number of different w(z), we find 
that the family of [TV, M] reconstructions which satisfy 

CIC < CIC min + k (4) 

where s is in the region of 0.2 and k = 5 yields very 
solid results. Typically the BIC produces models with 
few parameters and small errors and generally disfavours 
large variations in w(z) unless strongly warranted by the 
data; the AIC on the other hand renders more featured 
reconstructions, at the expense of larger errors. We find 
s = 0.2 works well for the models we present below, bal- 
ancing AIC and BIC. For alternative data sets, s and the 
choice of ft = 5 can be adjusted (e.g., n = 10 is more 
robust). 

The reconstructed cosmological constant discussed 
above is shown in Fig [2] In this case where s = 0.2, 
the reconstruction is good for z < 1, but degrades above 
z ~ 1.4. For s > 0.8, however, the fits at high z im- 
prove dramatically with the reduced freedom in the fit- 
ting functions; this is comparable in complexity to fitting 
a constant w, and so the errors are very tight. This result 
improves significantly on previous non-parametric recon- 
structions of A using SNAP-like data, such as in [TO] . 

Results In Fig. [3] we show the method in action for 
two very different types of w. One is a standard slow 
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FIG. 4: Reconstruction of w using the constitution SNIa. 



evolution, given by w = ^ { — 1 + tanh [3 (z — ^)] }; the 
other mimics the effective w(z) one finds in best fit 
void models of dark energy [19 , which we model by 
w = 0.2 + 1.8 exp [-(z-0.4) 2 /0.15j. We use s = 0.2 
and £ fiducial = £tme can gee reconstruction 
is impressive, with errors ~ 0.1. Above z ~ 1 the errors 
grow uncontrollably, leading to weak constraints despite 
the large number of SNIa in this range. This is because 
constant errors in \i lead to strongly divergent errors in 
w in a non-parametric setting. We find that similar fits 
for other w(z) models which have features on the same 
scales in z-space. We find that if the underlying w(z) 
contains sharp features such that its derivative is large, 
then a choice of s = 0.2 isn't sufficient, and the CIC 
needs to be weighted more heavily towards AIC (smaller 
s), which increases the errors. Choosing s = 0.05 with 
iV max = 10 lets us reconstruct a step-like w with a step 
of width ~ 0.1, or a Gaussian bump of about twice that 
width. Thus, we see that s (combined with the choice of 
largest A/", and a choice of k) sets the resolving scale for 
the reconstruction and must be treated as a prior, rep- 
resenting one's intuition of the complexity of w(z). For 
example, we can reconstruct w = — 1 to high accuracy by 
using 5 = 1, shown in Fig. [2] (right). This gives errors on 
w of less than 5% to z ~ 1. This situation is analogous 
to the choice of eigenfunctions which minimise the risk 
around A in [3]. 

As a final example we consider the w(z) we obtain 
from the Constitution data set [20 . To do this, we 
first find the best fit LCDM model, which gives ( = 
{0.32, —0.09, .653}, and we use this for both the resid- 
ual \i and for the model reconstruction. In Fig. [4] we 
show the constraints on w(z) with different choices of s, 
and using k = 5. This serves as an illustration of the 
method and the effect s has, but errors on £ have not 
been folded in for clarity. Using 8 = 1 yields tight con- 
straints on a par with those found assuming a constant 



w [20]. This results from using BIC as a selection cri- 
terion which yields conservative reconstructions, at the 
risk of smoothing over real features in the data and un- 
derestimating the errors. For realistic constraints from 
a given data set, simulated data of the same standard 
should be examined first to quantify sensible choices for 
s and n. 

Discussion We have presented a novel method to per- 
form a direct reconstruction of w(z), shown to be capable 
of constraining w to ~ 0.1 — 0.2 for z < 1 using SNAP- 
quality data. The method capitalizes on the use of prin- 
cipal component analysis to find orthogonal bases with 
which to fit the data, and uses a combination of infor- 
mation criteria to construct a family of reasonably good 
fitting functions. Then, by combining the different fit- 
ting functions we find the 'real' structure present in the 
data, and smooth over noise-induced features. All fitting 
is linear, so calculations and errors estimations are easy, 
with no exploration of parameter space required. Future 
work will incorporate errors from other parameters, and 
explore an iterative approach to finding /i res id based on 
the best fit obtained. 
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